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We present an improved version of the approximate scheme for generating inspirals of test-bodies 
into a Kerr black hole recently developed by Glampedakis, Hughes and Kennefick. Their original 
"hybrid" scheme was based on combining exact relativistic expressions for the evolution of the 
orbital elements (the semi-latus rectum p and eccentricity e) with approximate, weak-field, formula 
for the energy and angular momentum fluxes, amended by the assumption of constant inclination 
angle l during the inspiral. Despite the fact that the resulting inspirals were overall well-behaved, 
certain pathologies remained for orbits in the strong-field regime and for orbits which are nearly 
circular and/or nearly polar. In this paper we eliminate these problems by incorporating an array 
of improvements in the approximate fluxes. Firstly, we add certain corrections which ensure the 
correct behaviour of the fluxes in the limit of vanishing eccentricity and/or 90° inclination. Secondly, 
we use higher order post-Newtonian formulae, adapted for generic orbits. Thirdly, we drop the 
assumption of constant inclination. Instead, we first evolve the Carter constant by means of an 
approximate post-Newtonian expression and subsequently extract the evolution of l. Finally, we 
improve the evolution of circular orbits by using fits to the angular momentum and inclination 
evolution determined by Teukolsky-based calculations. As an application of our improved scheme 
we provide a sample of generic Kerr inspirals which we expect to be the most accurate to date, 
and for the specific case of nearly circular orbits we locate the critical radius where orbits begin to 
decircularise under radiation reaction. These easy-to-generate inspirals should become a useful tool 
for exploring LISA data analysis issues and may ultimately play a role in the detection of inspiral 
signals in the LISA data. 



I. INTRODUCTION 



Among the most promising and rewarding sources of gravitational radiation for the future LISA space-observatory 
are the so-called extreme mass ratio inspirals (EMRIs). These binary systems are formed as the result of the 
capture of compact stellar remnants by supermassive black holes in the nuclei of galaxies ,2| . As a result of two-body 
scattering in the surrounding cusp stellar population, such objects can end up on orbits that pass close to the central 
black hole. Occasionally these objects are captured by the massive black hole's gravitational "pit" with initial orbital 
parameters which ensure that the subsequent dynamics of the newly formed binary are almost entirely governed 
by gravitational radiation. The expected event rate for this scenario is quite promising (taking into account the 
capabilities of LISA), despite the uncertainties 0, 01 ■ 

The detection of EMRIs and the subsequent extraction for the system's parameters will rely heavily on the technique 
of matched filtering which requires an accurate theoretical model of the true gravitational wave signal. The 
computation of EMRI waveforms is an ongoing effort, which is currently limited by our inability to accurately describe 
the orbital motion of the small body in the Kerr spacetime under the influence of radiation reaction. The traditional 
way of describing the system is via the celebrated Teukolsky perturbation formalism ^ (where the system is modelled 
as a test-body in the field of a Kerr black hole) combined with the assumption of adiabatic evolution (see [3j0 
for recent reviews and further references). This latter property is well justified as a first order approximation: the 
system's extreme mass ratio /i/M ~ 10^^ guarantees that any influence of gravitational back- reaction will become 
significant at timescales much longer than any orbital timescale. In the adiabatic approach one first assumes that 
the small body is moving on a geodesic and computes averaged fluxes of energy E and angular momentum Lz for 
this orbit. These fluxes are used to update the parameters of the geodesic and the procedure is then repeated. This 
procedure works very well in special cases: orbits that are either equatorial, or inclined and circular. For those two 
families the change in the third integral of motion, the Carter constant Q, is either trivial or can be directly inferred 
from the other two fluxes. This is no longer possible for generic (i.e., eccentric and inclined) orbits. Dealing with 
them requires input from the more advanced framework of gravitational self-force computations (see |^ , |lClj| for recent 
reviews in the field). 

Generating an inspiral and the associated waveform is a computationally expensive exercise even at the level of 
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the Teukolsky formalism, let alone a self-force based calculation. Presently there are Teukolsky-based results for 
equatorial or inclined-circular orbits that show how a given orbit would evolve under radiation reaction , j 
and recently some initial results have become available for generic orbits jl^l- However, to date the only available 
full inspiral computation is the one by Hughes 'iS] for circular-inclined orbits. A computation along the same lines 
can be carried out for equatorial-eccentric orbits and full generic inspirals should appear in the near future. However, 
these will be computationally very expensive. 

Since present perturbative methods are either computationally demanding or incomplete it is highly desirable to 
have at hand a "quick and dirty" scheme for generating Kerr inspirals and waveforms, which can then be utilised in 
crucial data analysis computations for LISA 0. Glampedakis, Hughes and Kennefick jl^ (hereafter GHK) proposed 
a hybrid scheme for computing approximate generic EMRI trajectories. This hybrid scheme works by combining 
post-Newtonian radiation reaction formulae with a strong- field definition of the orbital parameters. It is able to 
reproduce many of the features expected from true inspirals and compares well to existing Teukolsky-based results 

[ni, [13, [la. 

However, there are some significant problems with the approach in certain regions of the orbital phase space. We 
describe in the present paper a set of fixes and improvements to the GHK formalism that address these problems. The 
most crucial improvement comes from the study of the limiting cases of nearly circular and/or nearly polar orbits. In 
order to guarantee a smooth orbital evolution, the exact fluxes need to satisfy certain consistency relations. By adding 
appropriate correction terms we can enforce the fluxes to satisfy these relations, thereby eliminating the pathologies 
of the GHK scheme in those regions of phase space. 

We also make improvements to the weak- flux expressions for E and L^. An immediate improvement comes from 
simply using available higher post-Newtonian ( PN) expressions. By combining Tagoshi's 2.5PN fluxes for equa- 
torial small-eccentricity orbits, Shibata et aZ.'s [T^l 2PN fluxes for circular small- inclination orbits and Ryan's [I3 
fluxes (which are fully accurate in the leading Newtonian and 1.5PN pieces ), we construct approximate 2PN fluxes 
which perform well for practically all eccentricities and inclinations and improve the strong-field behaviour of the 
GHK inspirals. 

A further modification is to change the prescription for the evolution of the orbital inclination, or equivalently 
for the Carter constant flux. GHK adopted the simple (but fairly accurate) rule of fixed inclination during inspiral. 
Here we improve on this by adding the next order spin-dependent correction, and allowing inclination to evolve. A 
flnal improvement to the fluxes comes by fltting functions to the results of Teukolsky-based computations for circular- 
inclined orbits ^12] . The combination of the above ingredients ensures the new hybrid scheme is applicable throughout 
parameter space, is qualitatively correct everywhere, and exhibits very good performance when compared to accurate, 
Teukolsky-based inspirals. 

The paper is organised as follows. In section ^] we give a brief overview of the GHK approach to computing 
inspirals. In sections IIIIIIIVI we describe two corrections that are required to give physically reasonable behaviour in 
nearly circular and nearly polar inspirals. These corrections must be included when any approximation is used for the 
energy and angular momentum fluxes. Then, in sections IVl IVII we provide new flux approximations that are accurate 
to higher post-Newtonian order than the leading order fluxes used in 16]. In section IVIll we illustrate the use of these 
expressions by investigating the properties of some example inspirals. Finally, in section rVIIII we briefly discuss the 
conservative component of the self- force which is not otherwise considered in this analysis. In the standard manner 
we adopt geometrised units G = c = 1. 



II. SUMMARY OF GHK HYBRID SCHEME 



In the GHK paper 0, the orbit is reparameterized in terms of an eccentricity, e, semi-latus rectum, p, and 
inclination angle, i. The parameters e and p are defined in terms of the turning points of the radial motion, which 
are determined by the roots of the radial potential 



R{r) = [E{r^ +a'^)~a L,] ^ - [r^ ~ 2 M r + a^) fi^ + {L^ ~ a Ef + Q 



(1) 



We use M and to denote the masses of the primary and secondary respectively. For a given value of the energy, 
E, z-component of the angular momentum, L^, and Carter constant, Q, the motion oscillates between a periapse, r^, 
and an apoapse, Tq. Once these have been determined from (Q, p and e are defined by 
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The orbital inclination is defined in terms of the Carter constant by 

Q ^ LI tan^ i. 



(2) 
(3) 
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In the limit a ^ 0, Q ^ + Ly = L"^ — L^, and definition ^ agrees with the usual notion of inclination. 

GHK constructed inspirals by evolving E, and Q. The energy and angular momentum were evolved using lowest 
order post-Newtonian fluxes, modified from Ryan jl^l 
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Ryan also provides an expression for the evolution of Q, but this was found to overpredict the change in inclination 
angle, when compared to accurate Teukolsky-based calculations for circular-inclined orbits fl^. For this reason, GHK 
used a "constant inclination" approximation^ to evolve Q 
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The energy and angular momentum fluxes (^HSJ are evaluated using the p and e given by definition This was 
found to work much better in the strong-field than re-writing the right hand sides of 101)-® using the Keplerian 
relation between E, and p, e and evaluating the resulting expressions for the orbit with the same E, and Q. 



III. CORRECTING THE BEHAVIOUR OF NEARLY CIRCULAR INSPIRALS 

The original GHK inspirals exhibit some unphysical features for nearly circular (e « 0) and nearly polar (l k, 1^/2) 
inspirals. These problems arise because the GHK approach uses post-Newtonian expressions to evolve E, Lz and Q, 
which are only valid to certain orders in the spin, eccentricity etc., but then computes the evolution of the eccentricity 
and semi-latus rectum using the exact transformation law, accurate to arbitrary order. As a consequence, sensitive 
cancellations which are needed to ensure reasonable behaviour in certain limits do not occur and this leads to some 
bizarre behaviour. This can be ameliorated somewhat by extending the flux expressions to higher post-Newtonian 
order (see section^), but the problems are not eliminated entirely. However, it is possible to make corrections to the 
fluxes that force the necessary cancellations to occur. We describe these corrections in this and the following section. 

In Figure ^ we show a selection of inspirals in the equatorial plane of a black hole with spin parameter a — 0.9M. 
Three inspirals are shown, with initial semi-latus rectums and eccentricities given by the pairs (20 M, 0.3), (15 M, 0.001) 
and (10 M, 0.001). The inspirals have been computed using the original GHK fluxes Q-® . The unphysical behaviour 
near e = manifests itself in two ways. For inspirals that start with comparatively large semi-latus rectum, but very 
small eccentricity, the inspiral moves rapidly away from circularity - in fact, the trajectories are almost vertical near 
to the e = axis. For inspirals that start with high eccentricity, the turnover point at which the eccentricity starts 
to increase occurs at a higher value oip/M than one would expect from Teukolsky-based calculations 0. There is 
also a strange tendency for the curves to attract - for these quite diverse choices of initial parameters, every inspiral 
plunges at nearly the same point in phase space. Changing from the original GHK fluxes to the higher order fluxes 
described in section alleviates the rapid increase in eccentricity prior to plunge, but the nearly circular orbits are 
still pushed away from circularity and the trajectories still converge at plunge. 



A. The cause of the problem 

The origin of this problem can be understood by considering the behaviour of the fluxes near e = 0. The rate of 
change of the eccentricity is given by the fluxes of energy, angular momentum and the rate of change of the inclination 



^ This rule is in fact exact in a spherically symmetric gravitational field, see Il6l for a discussion. 
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FIG. 1: Illustration of the problem with the hybrid fluxes for nearly circular orbits. We show three inspirals in the equatorial 
plane of a black hole with spin a — 0.9A/. 



angle, i (or equivalently, the flux of the Carter constant Q) by 
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The partial derivatives may be derived from the definitions of the radial potential and eccentricity (O 

R[r) = {E"^ - fi^) [r - ra)ir - rp){r - r3){r - ri) 
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Since the eccentricity enters this equation only as e^, it follows that E, and the fluxes E and must be functions of 
p, e^ and t. The post-Newtonian fluxes Q-© satisfy this condition. However, it also means that may be written as 
a function of E and L^. We anticipate that d{e'^)/dE is well behaved for all eccentricities, but de/dE = {d{e'^)/dE)/e 
will diverge like 1/e. Differentiation of the condition ^ yields the necessary partial derivatives 
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D{r) = -2 - E^)r^ + 3^1^ Mr^ - (a^ (^^ e^) + sec^ l) r 

+Llia.T? lM +{L^~ aEf M (10) 

The substitution e — e leaves the right hand side of these expressions unchanged, so as expected the derivatives 
are functions of only. The function D{r) is given by the derivative of the radial potential ^ dR/dr and so we can 
re- write it 



D{r) = -- (/i^ - E^) [- (r - r^) (r - rg) (r - + (r^ - r) (r - rg) (r - r^) 
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In the limit of small eccentricity at fixed p and t, we deduce 
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where Ni{p, l) = Ni{ra) = N2{rp) for 6 = and similarly for A^2, ^3- 

It is now clear from (jS)) that e cx 1/e for small e, unless the fluxes satisfy a consistency condition in the limit e — > 

A^i(p, £; (p, 0) + N2ip, l) L, (p, l, 0) + Nsip, i) L {p, i, 0) = (13) 

We will see in section Hvl that it can be better to evolve the Carter constant, Q, instead of i. Using the definition Q 
to write i in terms of Q and , (|13|l becomes 

Ni{p, i) E {p, i, 0) + N4p, i) {p, i, 0) + N^ip, l) Q (p, 6, 0) = (14) 

where 

N4,{p,i,) ^ {2Mp-p^)L,-2MaEp, N^ij), l) ^ {2 M p - p^ - a^) /2. (15) 

It is understood that the A^^'s are evaluated for the circular orbit with the specified p and l. The GHK fluxes do 
not satisfy this relation, which explains the bizarre behaviour seen for small eccentricity. 

What does the condition (|13|l or (I14|) mean physically? For simplicity, we first examine equatorial orbits before 
extending to inclined orbits. In the equatorial plane, i = by symmetry, and this condition is satisfied by the GHK 
scheme In the circular equatorial limit, the energy and angular momentum are given in terms of the semi-latus 
rectum by 

E _ 1-2 {M/p) ± (a/M) {M/p)i ^ ± (P\ ^ 1 T 2 (a/Af) (M/p)^ + {a/Mf {M/pf ^^^^ 



Y^l - 3 {M/p) ± 2 {a/M) {M/pY ' ^1-3 {M/p) ± 2 {a/M) {M/p) 

Where the upper (lower) sign is for prograde (retrograde) orbits. We find that condition (|13|l reduces to 

E{p, O/tt, 0) = ± 3 4 (p, 0/^, 0) = n^{p) L,{p, O/tt, 0) (17) 

where ri<^(p) is the azimuthal velocity, d(j)/dt, of the circular orbit. Expression (|17|l is just the condition that circular 
orbits remain circular under radiation reaction, which must be true for the actual gravitational waves emitted by the 
orbit The general condition (|13|l is similarly just another way to write the circular goes to circular rule |2lj . By 
forcing the leading term to cancel from our e expression, we ensure that e cx e and therefore that e = for circular 
orbits. 



6 



We can compute the difference E — Vl^p for circular equatorial orbits under the GHK scheme. Expanding in 
powers of the spin parameter, a, we find 
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The cancellation occurs correctly to linear order in the spin. This is expected since the fluxes were derived in a 
low-spin limit. If we use the higher order fluxes (|37ll - (|38|l described in section|3 we similarly find 
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In each case, the cancellation occurs to the post-Newtonian order at which the fluxes were derived. However, when 
we use the fluxes to evolve inspirals, this is no longer accurate enough since we are using definitions of eccentricity 
and semi-latus rectum that are accurate to arbitrary order. Since we know that the circular goes to circular rule must 
hold for real inspirals, we can improve our approximate scheme by modifying our fluxes to enforce this condition. 



B. Correcting the hybrid fluxes 



The condition (|13(l is a relationship between the fluxes for circular orbits. The easiest way to impose it is to specify 
two of the flux expressions and use H13|) to determine the other. Since this relation applies only to circular orbits, 
we adopt the approach of leaving the eccentricity corrections largely unchanged, and correct only the circular piece 
of the fluxes. However, we must include the usual (1 — e^)^ factor in front of the near-circular correction to ensure 
reasonable behaviour at high eccentricity. As the eccentricity approaches the parabolic limit, e = 1, at fixed semi-latus 
rectum, we expect the total change in E^ L^, Q and i to be approximately constant, since the strong-field part of 
the orbit^ where most of the radiation is lost, does not change significantly. However, the orbital period diverges like 
(1 — e2)~2 ^ and therefore the orbital averaged flux must tend to zero in the same way. This is the origin of the (1 — e^)^ 
prefactor in Q-©, and we must also include it in the circular correction. We cannot reliably say anything about 
further eccentricity corrections, since such corrections are higher order than the available post-Newtonian results. Our 
approach ensures that the eccentricity pieces are asymptotically correct and the new fluxes have the right leading 
order behaviour at high eccentricity. It may be possible to get more accurate expressions by applying the correction 
to the other eccentricity terms as well, but a priori we have no justification for doing this. 

If we chose to fix E and Lz by expressions and use H13|) to determine i, we would find that Z ^ in the 

equatorial plane, which is unphysical. For simplicity, we want to use the same correction throughout parameter space, 
so we will modify either E or ■ In Figure |2 we show new equatorial inspirals computed by imposing the fix in 
two different ways. For the solid lines we use the post-Newtonian E expression Q), and correct the circular piece of 
the post-Newtonian using H13I) . For the dashed lines, we leave the post-Newtonian expression ((SJ unchanged, 
and correct the circular piece of the E expression. We note first that it does not make a great deal of difference 
which way the fix is implemented. The trajectories differ a little near to plunge, but are largely the same. The most 
striking feature of Figure |21 is how much better the inspirals now look. The orbits are no longer pushed away from 
circularity and the turnover in the eccentricity occurs much later - both features that are expected in true inspirals. 
The different trajectories also plunge at somewhat more different points. There is some excess growth in eccentricity 
close to plunge, but this is due to the failure of the GHK flux expressions Q-®. Using higher order expressions (see 
sectional, eliminates this feature. We conclude that it does not make too much difference how the fix is imposed, but 
imposing the fix gives a significant improvement in the qualitative behaviour of the inspirals deep in the strong-field 
region of the spacetime. 

In the limit of a nearly polar orbit, we find that eidej dL^\E,L and e{deldi)\E,L both diverge like l/Lz, while 
e (9e/9i?)|L.t is finite (this will be discussed in more detail in the next section). As a result, if we specify Lz and 
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FIG. 2: Improved inspiral trajectories after implementation of the fix described in the text. We show inspirals in the equatorial 
plane of a black hole with spin a = 0.9M. The solid lines correspond to using the post-Newtonian E expression Q to correct 
the angular momentum flux, L^- The dashed lines correspond to using the post-Newtonian Lz expression @ to correct the 
energy flux, E. 



i and use H13|l to correct the circular piece of E, the correction diverges in the limit l tt/2. Therefore, when 
evolving l one should use E and i to correct L^. However, if we evolve Q instead of t, the functions e {de/dLz)\E,Q 
and e {de/dQ)\E,L are both finite at the pole and there is no problem. However, e {de/dLz)\E,Q vanishes for certain 
retrograde orbits. Using E and Q to correct would therefore lead to a divergence at this point (unless E and Q are 
constrained to cancel appropriately). Thus, when evolving Q, it is better to use Lz and Q to correct E. The modified 
E is 
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In this, {E,Lz, Q)ghk refer to the prescription of the energy and angular momentum fluxes being used, e.g., (O, 
in this case. If the fluxes are modified (e.g., to the 2PN expressions given in sectional the correction is applied as 
above, with {E)chk replaced by {E)2PN etc. 



IV. CORRECTING THE BEHAVIOUR OF NEARLY POLAR INSPIRALS 



The GHK hybrid scheme also exhibits some pathological behaviour for orbits that are nearly polar, i.e., with 
t « 7r/2. In Figure 121 we show inspirals for four different initial inclination angles, computed under the original GHK 
scheme. The inspirals have two strange properties. Firstly, the transition passing across the pole is not smooth - 
there is a discontinuity between trajectories with l slightly less than 90°, and those with t slightly above 90°. In fact, 
it is not possible to compute trajectories in this approximation for i slightly less than 90°, since it predicts an increase 
in p. Secondly, orbits with high prograde inclinations rapidly circularise, which is unphysical. This latter effect 
looks like a manifestation of the problem e oc 1/e for small e discussed in the previous section, but the coefficient of 
proportionality is now positive. We might hope that the correction for near-circular orbits discussed above will fix this 
problem. Figure^shows the same set of inspirals, computed including the near-circular correction l|20|) . As expected, 
the correction does fix the rapid circularization problem, but the transition across the pole is still discontinuous. 
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FIG. 3: Inspirals into a black hole of spin a = 0.9 M for initial eccentricity e = 0.3 and initial semi-latus rectum p = 10 M. 
These inspirals were computed using the fluxes Q-©, but without the small-eccentricity correction. The inspirals are shown 
for initial inclinations of t = 75° , 80° , 85° , 100° and 105° . 



The remaming problem near l — tt/2 can again be understood by considering the behaviour of e (or p) in that 
limit. The derivatives de/dL^ and de/di are both singular, diverging like sect as t — > tt/2. This is clear from looking 
at N2{r) and N^{r) in that limit (|10|l . In fact, in terms of the Carter constant Q of the polar orbit, we find 



N2{r) ~ y/Q{~r^ + 2Mr - a^) sect, 7V3(r) « Q (-r^ + 2 Mr + a^) sect 



(21) 



The divergence of e, which corresponds to a divergence in fp and is unphysical. It will be avoided if the singularities 
in the and t fluxes cancel appropriately, which implies 



e, t : 



(22) 



There is nothing special about polar orbits which should allow us to specify an additional condition, so where does 
(|22|) come from? This is answered by looking at the derivative of the Carter constant 



Q = 2Lz tan lL^ + 2L^ tan t sec 1 1 = 2 ^JQ 
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The condition H22I) is equivalent to requiring that Q is finite in the limit t — > 7r/2. The problem has arisen because 
we have chosen to evolve the inclination angle rather than evolve Q directly. A real radiation field will generate a 
finite change in the Carter constant, and the rate of change of t computed from this will automatically satisfy H22|l . 
The main thing this tells us is that the constant inclination approximation, t = 0, is unphysical - polar orbits lose 
angular momentum, which drives them to become retrograde. An alternative prescription for the evolution of t is 
therefore required. If we specify and t at some post-Newtonian order, the relation H22I) will not hold exactly, as 
^/Q is accurate to arbitrary PN order. We can account for this by adding a correction to Lz, as for the near circular 
correction. However, since this correction is only required at the pole, it is difficult to extend it smoothly to non-polar 
orbits. The need for a correction can be eliminated entirely if we evolve Q instead of t, i.e., we choose to specify Q 
and derive t from this. 

Ryan's expression [l9l | for Q, written in terms of our orbital elements is 
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FIG. 4: The same inspirals as in Figure ^ but now computed using the fluxes (0J-||SJl and the small-eccentricity correction 
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We can use this in conjunction with the fluxes to compute inspirals. Figure \^ shows the set of inspirals 

depicted in Figures |31 and ^ but now evolved using (|24|l instead of the constant inclination approximation. The 
resulting inspirals now change smoothly as the inclination angle is increased. However, the change in t over the course 
of the inspiral is much too large when compared to the results of Teukolsky computations for circular orbits |15| . This 
is the reason GHK chose to evolve t rather than Q. In section IVll wc will show how to improve the approximation for 
Q to give more reasonable evolutions. 



V. INCLUSION OF HIGHER ORDER POST-NEWTONIAN FLUXES 



The corrections discussed in the previous two sections must be applied in any version of the hybrid scheme, and 
ensure that the inspiral properties are physically reasonable. However, the expressions used to evolve E, and Q 
can be improved, to give better agreement with perturbative calculations |12. .13^ .15j . 



A. Equatorial orbits 

A natural way to extend GHK is to replace with higher order expressions for the evolution of E and 

Lz- Tagoshi [13 derived 2.5 Post-Newtonian (FN) fluxes for the case of equatorial and eccentric orbits under the 
assumption of small eccentricity (Tagoshi kept up to C(e^) terms). The assumption of small e is not as restrictive as 
it initially sounds, as can be demonstrated by comparing inspirals generated using the GHK fluxes I©-® to those 
generated using GHK truncated at C(e^). This comparison is shown in Figure Eland indicates that inspirals can be 
faithfully reproduced using the e-truncated fluxes provided e ^ 0.8. We expect the same to be true for the higher 
order FN fluxes that we construct below. 

In his calculation, Tagoshi |^ describes the test-body's orbital motion with the parameters ro, ej and v"^ — M/tq. 
The parameter tq is the radius at which the potential, Vr, governing the radial motion has a minimum. Then, ej is 
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defined by the requirement that r = ro(l + ej) is a turning point for the radial motion, i.e. Vr{r = ro(l + ej)) = 0. It 
turns out that the radial motion r{t), to C(e^) accuracy, is given by, 



where 



r{t) - ro[l + e,A'\t) + ey^\t) + 0{e% 



A^\t) = (73(1 - cos ri^t) + (74(1 - COS 217^^). 



(26) 

(27) 
(28) 



The quantity Qr is the angular frequency associated with the radial motion, while 93,94 are functions of vq and the 
hole's spin q — a/M (we write explicitly only the former). 
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From equation H26(l we can immediately see that the apoastron and periastron are given by, 

ra = ro(l + ej)+0(e3), 

rp = ro{l-ej+2q3e^j) + 0{e^j). 



(29) 
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The averaged fluxes at infinity are given by [l7j | , 
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FIG. 6: Approximate equatorial inspirals for an orbit with initial parameters p = 20M, e = 0.8. The dashed curve was produced 
by use of the fuU-e Ryan fluxes while for the solid curve only 0{e^) terms were retained. The black hole spin is a = 0.5M. 
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In order to adapt these fluxes to the GHK scheme we first need to rewrite them in terms of our parameters p and 
e. Assuming a small eccentricity, e <C 1, we have, 
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Combining these with equations (|30|I - H31() leads to 



e + e^q^ip), 



1 = l[l + e\q,ip)-l)]. 



(35) 
(36) 



It is now straightforward now to rewrite the fluxes H32|l in terms of p and e. To enhance accuracy, we adopt the 
approach 'where possible, include higher order terms in e' (for example, by comparing with Ryan's expressions). In 
particular, we must have the factor (1 — e^)^ to ensure the behaviour is qualitatively correct for high eccentricity, as 
discussed earlier in the context of the near-circular correction. We find. 
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where the various e-dependent coefficients are, 



51(e) 
53(e) 
55(e) 
57(e) 



73 

^+24' 

1247 
"336" ^ 

44711 
9072 ^ 

8191 



2 37 
+ 96^ 

9181 n 



672 

172157 
^ 2592 ^ 

44531 o 



672 



336 



59(e) = 1 + -e^ 



5ii(ej 
513(e) 
515(e) 



1247 425 



336 

44711 
9072 

8191 



336 ' 

302893 
^ 6048 ^ 

48361 o 



672 



1344 



510(e) -( — I 511(e) +71- ( — I 312(6) 



5/2 



516(e) 



52(e) 
54(e) = 4 



73 
12 

1375 



823 



949 



58(e) 



514(e) 
516(e) 



48 


-e , 


33 


359 2 




e 


= 16 + 


32 ' 


3749 


5143 




e 


336 


168 


61 


119 2 


= 12 + 


+- 




97 2 


= 4 + 




33 


95 2 






~ 16 ' 


^16^' 


417 


37241 


~ 'W 


672 



183 

1 

32 



— 513(e) 
p J 



(38) 



491 
192^ 



(39) 



It is a well known fact that PN expansions are characterised by poor convergence, that is, a higher PN order does 
not necessarily mean more accurate result. The same behaviour is found in these fluxes too. After some testing we 
found that the optimal order is 2PN (in agreement with the literature on the subject). This is illustrated in Figure [7| 
where we compare 2PN with 2.5PN fluxes. Note that in this and subsequent figures, we are including the near-circular 
correction (|20|l and, for inclined inspirals, are evolving the inclination by prescribing Q. 

Figure IHl illustrates a set of inspirals computed using the revised fluxes (truncated at 2PN level). These 

represent equatorial inspirals into a black hole of spin a = 0.8M. This is one case for which the original GHK scheme 
did not perform accurately. The revised inspirals are a significant improvement over the GHK results, particularly 
near the separatrix. Comparison of the numerical values of and E with the results of Teukolsky-based calculations 
[13| also indicate we are making a significant improvement. 



B. Extension to non-equatorial inspirals 



The results presented in the previous section apply only to equatorial orbits, but ultimately the goal is to improve 
the inspiral of generic (inclined and eccentric) orbits. The only available higher order PN expansion for inclined orbits 
has been derived by Shibata et al. _,18j for the case of circular orbits with small inclination. Shibata et al. use the 
inclination parameter 

and assumed that this is a small number. The relation between y and our inclination angle l is, 

'°'''=l + y[l+a2(l-£;2)/L2]' (41) 
which for small t, y and to 2PN accuracy becomes, 

V^i\\- a\\ - i?2)/L2] ^ ,2[i _ {alMf{Mlvn (42) 



13 




I I I I I I I I I 

4 6 8 10 p^i^ 12 14 16 18 20 

FIG. 7: Determining which PN order is more rehable: the figure shows inspirals with initial parameters p = 20Af, e = 0.5 and 
spin a — 0.5M generated using the fiuxes (1371 . truncated at 2.5PN (dashed curve) and 2PN (solid curve). It is clear that 
2PN is the best choice. 



The fluxes derived by Shibata et. al. are given by. 
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At 2PN accuracy, we find that 1 — y/2 = 1 — t^/2 + / 2{a / {M / pY . Term by term comparison (where possible) 
between expressions and Ryan's fluxes 0)-® suggests the correspondence 1 — y/2 cosi and 1 — 3y/2 — > 
— 1/2 + 3/2 cos^ L when one goes from small to arbitrary inclination for these terms. For the 2PN terms, we expect 
E to contain terms proportional to 1 and cos^ t, while should contain pieces proportional to cost and cos'^ i [2^ . 
This is born out by fits to circular Teukolsky computations {^^^ section IVI C 1|) . 

The cosine terms contribute in the equatorial plane and therefore, using (|37|l - l|38|) . we know the e-dependent factors 
that multiply them. However, this is not true for the sine terms. Consequently, we expect terms that contain sin b to 
be incomplete and of modest accuracy. Nevertheless we have included them in our fluxes (and verified that for orbits 
not close to the equator these terms play an important role). Putting together all of the above, we extend the fluxes 
l|57jl - l|55|l to generic orbits as follows 
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As discussed earlier, we have included higher order eccentricity terms where possible and truncated the series at 2PN 
level. We must also remember that we do not actually use the circular pieces of the E expression (|44|l . but determine 
these using the circular fix H2U|) . 

Some numerical inspirals produced by incorporating these fluxes in the GHK scheme are shown in Figure |51 along 
with the original GHK inspirals for comparison. These are inspirals into a black hole with spin a = 0.9M for two 
different inclinations - l — 30° and l — 60°. We are evolving Q using Ryan's expression H24|l . As for the case of 
equatorial orbits, we find that the new fluxes produce inspirals that are more physically reasonable than the original 
GHK results. 



VI. AN IMPROVED PRESCRIPTION FOR THE CARTER CONSTANT FLUX 



The constant inclination approximation ((TJ is a "spherical" approximation in the sense that it is exact in the limit 
a 0. In section HVI we saw that i = cannot be used to evolve nearly polar orbits. High inclination orbits offer 
the orbiting body an opportunity to more effectively probe the aspherical nature of the Kerr metric, thus it is not 
surprising that the "spherical" rule fails in such cases. Unfortunately, the available (and approximate) expressions for 
computing Q for generic orbits are extremely limited. Ryan's weak-field expression [T9l | gives rather poor results. As 
discussed in |0, the bulk of the change in Q is reliably captured by equation {T)). Hence, we will attempt to add a 
correction to this expression in order to improve the inspirals. 
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FIG. 9; Inspirals generated with the help of the new fluxes 14411 . I45I I (sohd hnes) compared to the original GHK inspirals 
(dashed lines). The initial parameters are p — 20M, e — 0.6, t — 30° and t = 60° for black hole spin a = 0.9M. 



A. What is the connection between the "spherical" Q and the Kennefick-Ori formula ? 

Kennefick and Ori give an expression for the evolution of the Carter constant in terms of the theta component 
of the self-force. Expressing the Carter constant as 

the Kennefick-Ori formula (equation (11) of j2]|) gives, 

Q^H^EE + H^L^L, + ^—Fe, (48) 
where S = 7'^ + cos^ 6. Writing this explicitly, 

Q = ~2a^Ecos^e E + 2L,^cot^e L, H —Fe. (49) 

This formula and the spherical lITIl have both been discussed and used in calculations in the extreme mass ratio 
literature over recent years [T^ ITo. l2ll| . However, as yet there has been no discussion on their relation. 

In the a — > limit, the first term of (|49|) vanishes, E ^ r^, the Carter constant becomes Q = L\^^ — (where 
Ltot is the total angular momentum) and the inclination angle i becomes the true inclination of the orbital plane, 
tan^ L = Q / L\. Taking the y-axis to lie in the orbital plane, without loss of generality, the Boyer-Lindquist coordinates 
of the particle at any point of the orbit obey the relation 

/(0, (f) = cos 9 — sin 6 cos tan l ~Q . (50) 

Symmetry ensures that the self-force has no component perpendicular to the orbital plane and therefore 

sec 6 tan ( 



u-n = = F- n ^ _ = _ = — Z. Z. = sin^ y/tan^ t - cot^ 9 . (51) 

In this, ria ~ df /dx"' is the normal to the orbital plane and — da;"/dT is the particle four-velocity. The third 
term in equation (|49() thus becomes 

2i:u'^ 2r^u^ F^ _ ( u'^\^ ^u't' F'l' (tan^ t - cot^ 6*) 



Po- - 2 r-—-=2^- '-u,F,. (52) 
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But, Ucf, = Lz and = dL^/dr, so F^/u* = and H49|l becomes 

Q = 2 (cot^ 6* + tan^ l - cot^ = 2 tan^ lL^L^ ^2Q^ , (53) 

which brings (|49(l to the form (0. Switching back on the spin makes (I49II depart from its spherical value and 
consequently l to change. It is easy to deduce that some leading order (with respect to the spin) terms will be 
provided by the third term in (|49|l . This is because the leading order change in Fg is linear in a, as can be found from 
Ryan's expressions for the self-force T^. The first term, on the other hand, contributes only at 0{a^) and beyond. 



B. An improved formula for Q 



If we expand the spherical formula iQ at the PN order of Ryan's expression, we find 
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where /4(e) = 61/24 + 63e^/8 + 95e''/64 as in equation ®, and q = a/AI as before. We see that the 0{a) piece of this 
expression is not the same as in equation H24(l which by construction includes all terms of this order. We conclude 
that there is an C(a) piece missing in the spherical Q formula. Written explicitly. 
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This piece represents the leading "aspherical" contribution to Q. This gives rise to an evolution in t. We note, 
however, that this term is divergent at the pole, since it is proportional to 1/ cost. This is a manifestation of the near 
polar problem discussed in section Hvl In reality, Q will be finite at the pole, so physically L must be what is required 
to cancel all of the divergent pieces of the "spherical" Q. We can therefore derive an improved prescription for Q by 
expanding the spherical formula iQ and removing all terms that diverge at the pole. In this way, we derive from the 
2PN angular momentum flux (|45|l a new expression for Q of the form 
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This method for prescribing Q removes the divergences at the pole, but is incomplete, since it tells us nothing about 
the pieces of L that vanish at the pole (or equivalently, the corrections to the non-divergent pieces of Q that come 
from i ^ 0). Nonetheless, this is an improved prescription for the evolution of Q. In the next section we derive a 
further improvement to the circular piece of the Q expression. 



C. Fitting to Teukolsky data 

1. Circular orbits 

The evolution of circular-inclined orbits has been computed accurately by solution of the Teukolsky equation 0,^3- 
We can use these results to improve the circular piece of our expression for Q. It turns out that the best results are 
obtained by fitting functions to the Teukolsky Lz and i data, and deriving Q from these, via expression H23|) . 

Using data provided by Scott Hughes, we find that a good fit to the angular momentum flux Lz for circular orbits 
is given by the function 
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where = -10.7420, 

= -33.7090, 
cl = 2.37258, 
cl = -7.61034, 
cl = 306.119, 
4 = 224.227, 
cg = -0.645000, 
/a" = 483.266, 
/I = 82.0780, 



cl = 128.878, 
cl = 40.9259, 
= -490.982, 

cl = -5.13592 
/2 = -1325.19, 
f?^ = 301.478, 



P J 
-9.07738, 

: 60.9607, 
-0.715392, 
cl = -475.465, 
c\ = -347.271, 
= -9.00634, 
Cg = 47.1982, 
= -219.224, 
fl = -904.161, 



: -1.42836 

40.9998, 
= 3.21593, 
cg = 12.2908, 
cg = 886.503, 
cl = 91.1767, 
= -283.955, 

■b 
3 



/g^ = 634.499, 



-271.966, 



= 10.7003, 

-0.348161, 
5.28888, 
cl = -113.125, 
c« = -25.4831, 
cl = -297.002, 
/i = 736.209, 
/4 = -25.8203, 
= 827.319. 



Similarly, a good fit to the evolution of t is given by 
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(58) 



The pieces of and L that arc non-vanishing at the pole were constrained to cancel so that the derived Q is finite 
there. The "d" coefficients describe these terms. The spin and inclination terms are polynomial in q and cos l and take 
the form of the next post-Newtonian terms that should enter the flux expressions. The j>dependence of the fit starts 
at 2.5PN order, but the factors of p are not consistent with the post-Newtonian orders of the spin and inclination 
terms they multiply. The "c" coefficients were obtained by a fit to data in the range p = 5M to p = 30M. A fit that 
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was consistent with post-Newtonian expansions was also derived, but this exhibited pathological behaviour when it 
was used at low p. The above form of the fit is virtually identical in the range over which the fit was derived, but 
does not show the same problems for p M . Data was also available for some orbits with p = 3M. The parts of 
the fits parameterised by the "f" coefficients were derived to reduce the errors in the first fit at p = 3 A/ and p = 5AI. 
The p-dependence was taken to be of a simple form that went to zero asymptotically at least as quickly as the rest 
of the fit. This ensured the quality of the fit was not degraded at large p. Note that several of the "f" terms in the 
Lz fit have the same form as the "c" coefficient terms. These are separated in H57|l only because they were derived in 
these two different ways. 

The fit matches the data to an accuracy of < 3% at all points. It should be trustworthy in the region p > 3M, but 
it will be less accurate as p -^M since data was not used in this range. However, the fit is such that the behaviour 
in this regime is qualitatively correct, namely < for near equatorial prograde orbits and i > for all orbits. 
Inspirals generated using this hybrid scheme will therefore evolve in a reasonable way, but the evolution will not be 
entirely accurate as the particle gets close to the central black hole. As this last stage of inspiral is very rapid anyway, 
the loss of accuracy in this regime should not affect the use of the hybrid scheme for exploring the inspiral problem. 
Moreover, the flux comparisons given in section FVII Bl indicate that the hybrid scheme does pretty well even for very 
strong-field orbits. 

As discussed in section Hvl we want to specify Q rather than L in the hybrid scheme. A suitable expression for Q 
is obtained by substituting expressions H57|) and (|58|l into equation H23|l . By including the factor sin^ i-j \fQ in the I 
fit, we ensure that the divergent terms cancel precisely and Q is finite everywhere. The procedure of fitting I and 
deriving Q may seem more convoluted than fitting for Q directly, but it ensures that the evolution of the inclination 
angle is always sensible, and thus generates more physically realistic inspirals. We have not presented a fit for E here, 
but this is not necessary since the evolution of energy for circular orbits is determined precisely by the evolution of 
Lz and Q when we apply the circular orbit correction described in section UTTl 

This procedure has allowed us to improve the hybrid fluxes for circular orbits, but we have not yet improved 
the eccentricity-dependent terms in the fluxes. However, as discussed earlier, we must include the usual (1 — e^)^!"^ 
prefactor to ensure reasonable behaviour in the limit e — s- 1. The resulting expression for Lz is 

[l\ = (l-e2)i [(l-e2)-i (4) (p, i, e, a) - (l,) (p, 0, a) + (l,) 1 (59) 

V /mod L V /2PN V /2PN V / fitJ 

in which the subscript "2PN" refers to expression and the subscript "flt" refers to expression H57f) . We obtain 
the expression for Q/^/Q in the same way 
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The final expression for E is still obtained from equation (|2()|l and expression H44|l 



2. Eccentric orbits 



The expressions quoted above are missing some eccentricity-dependent pieces. The omission of these eccentricity 
terms manifests itself as some unphysical behaviour for inspirals generated using the hybrid scheme. For prograde 
equatorial orbits with radius p = 2M around a black hole of spin a = 0.99M, we find that the magnitude of Lz 
becomes smaller as the orbital eccentricity increases, in stark contrast to the behaviour of Teukolsky generated fiuxes 
(see Table nj. For eccentricity e > 0.3 we find Lz > and p > which is unphysical. Glampedakis and Kennefick 
|13| tabulate Teukolsky-based fluxes for equatorial eccentric orbits and these can be used to derive an additional 
eccentricity-dependent correction to the fluxes. However, the inclination dependence of this correction is unknown 
since the data is prov ided for equatorial orbits only, and it is difficult to constrain the p dependence with the limited 
data provided in [l^j. Recently, Sago et al. published new 2PN results for the energy and angular momentum fluxes 
from orbits with small eccentricity and small inclination p^ . After rewriting the expressions in our coordinates, these 
results confirm the extrapolations that lead us to equations 1)44(1 - (|45|l and also give the eccentricity dependence of one 
of the 2PN terms that we could not otherwise derive. We deduce that equations ((44|I -H45 |I should only change in the 
final term, with the modifications 
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When these expressions (adding the consistency corrections and the fit for the circular flux) were tested by comparison 
to Teukolsky-based results, they were found to perform less well than the version of the hybrid scheme described above. 
It is not clear why this is the case. The new results in [2^ are inconsistent at 2.5PN order with the results in 
Although we do not use the 2.5PN order results, it is possible the inconsistency is indicative of some other error which 
might also have affected the new eccentricity terms above. However, it may just be that the new expressions do not 
perform well because we are using a weak-field expansion (the PN expansion) and applying it in the strong-field. The 
motivation for developing the hybrid scheme was to find a set of expressions that reproduce accurate, Teukolsky-based 
results as closely as possible. In that spirit, we recommend ignoring these new corrections and using the expressions 
quoted in sections I VI VI C l1 instead. Once the new results have been more thoroughly tested, and perhaps augmented 
with further fits to Teukolsky-based data, it might be possible to improve our expressions further, but this will be 
considered in future publications. 

Recent work by Sago et al. , based on a paper by Mino '25'| , provides an expression for the time evolution of the 
Carter constant, Q, which uses the same Teukolsky variables needed for the computation of the energy and angular 
momentum fluxes. The paper makes use of this expression to derive a 2PN formula for Q as well. The analogue 
of equation (|56|l is found to be 
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Note that Sago et al. actually give a 2PN expression for Q, not QI\fQ as given above. We quote the above result, 
obtained by expanding \fQ at 2PN order and comparing to Sago et al.'s results, for consistency with (|56|l . Also note 
that Sago et al. 's results were computed to linear order in i? , so it is not possible to determine any pieces proportional 
to sin^ L in the final term of equation H61|l . although such terms are present in expression H56|l . 

Equation l|61|) shows amazing agreement with our previous result, considering the way in which the latter was 
derived and this gives us some confidence in our approach. At present, there is no strong-field Teukolsky-based data 
for Q, so we have no way of determining if equation (|61|l performs better than (|56|l . Due to the inconsistencies in 
j2.ll | mentioned above, at present we suggest ignoring H61(l and using l|56|l (more precisely, (|56|l augmented with the 
circular fit H6U|l ). As Teukolsky data based on the new formula for Q j2j| becomes available, we will be able to assess 
the alternative expressions for Q and determine which gives the best results in the hybrid scheme. 

The results given in the following section are based on equations H44|l - H45|l and (|56|l . augmented with the near- 
circular correction H13|l and the fits to circular Teukolsky data H57|l - (j6()|l . We recommend using this form of the hybrid 
scheme in calculations, although one should bear in mind that the eccentricity pieces in these fluxes are incomplete 
and the fluxes are unphysical in a few small regions of parameter space. However, the unphysical behaviour appears 
only very close to plunge, thus the hybrid scheme in its current form should be adequate for most applications. Once 
Teukolsky-based data is available for a wider range of generic orbits it will be possible to compute a significantly more 
accurate fit to the missing eccentricity-dependent pieces. As we saw earlier, it is likely that fitting only the pieces 
of Teukolsky-based fluxes will be sufficient to give accurate inspirals and this will be pursued in the future. 



VII. INSPIRAL PROPERTIES 

A. Stability of nearly circular orbits 

Kennefick |23| examined the stability of nearly circular orbits in the equatorial plane using Teukolsky-based calcu- 
lations. This involved expanding the eccentricity derivative e near e = 0. As described above, we expect e w /(p) e 
near e = 0. For orbits very close to the separatrix, e > as e — > 0, while for orbits with larger periapse, e < in that 
limit. It is possible to find the point where the transition in behaviour occurs, i.e., where /(p) = 0. The locus of points 
where e = is important for determining the global properties of inspiral trajectories, and the root /(p) = is the 
point where this curve intersects the e = axis. Figure |51 indicates that the improvements to the hybrid scheme move 
the point where e changes sign in a given inspiral much closer to plunge when compared to the original GHK scheme. 
This is in keeping with the results presented in 20] . We can compute the root /(p) = as a function of the black 
hole spin, a/Af, for the new hybrid fiuxes and compare this to Kennefick's results (these were corrected in [T5 | and 
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FIG. 10: Location of the critical radius rcrit, at which the e = curve intersects the e = axis for equatorial orbits, as a 
function of spin. The pluses mark points computed using Teukolsky fluxes by Kennefick |20l| . while the crosses indicate points 
computed using the hybrid scheme described in this paper. The agreement between the two curves is within ~ 5% for all spins. 



we use the corrected values in this plot). This comparison is shown in Figure ITUl We find that the agreement in the 
location of the critical radius is extremely good, within ^ 5% for all values of the black hole spin. This is particularly 
remarkable since the location of the critical radius depends on the leading eccentricity terms in the fluxes, and we 
have already seen that the 2PN expressions used are somewhat inaccurate close to the separatrix. The structure of 
the inspiral phase space depends significantly on two curves - the separatrix at which orbits plunge, and the critical 
curve e = 0, where inspirals turn up. The hybrid scheme includes the correct separatrix by construction and Figure [TUI 
indicates that the new hybrid fluxes recover the critical curve very well, so we can have some confldence that inspirals 
generated under this improved hybrid scheme are qualitatively correct. 



B. Flux comparisons 



The hybrid flux expressions allow us to compute the rates of loss of energy and angular momentum from an 
arbitrary geodesic orbit. The true fluxes are given by solution of the Teukolsky equation, and we can compare the 
hybrid results to existing Teukolsky-based results available in the literature. This comparison is shown in Tabled 
The table is divided into three sections. The top section is for circular and inclined orbits, for which the Teukolsky 
data was provided by Scott Hughes IT3|. The middle section is for eccentric and equatorial orbits, for which the 
Teukolsky data was taken from Glampedakis and Kennefick The bottom section is for generic orbits, and the 
Teukolsky data was taken from Drasco and Hughes In each case we tabulate the hybrid and Teukolsky-based 
fluxes of energy and angular momentum, plus the rate of change of inclination angle. Drasco and Hughes use the 
"spherical" formula (TJ to evolve the inclination angle, which we have seen to be inaccurate. The Teukolsky results 
should not therefore be regarded as more accurate for that comparison. 

Overall, the agreement between the hybrid and Teukolsky results is remarkably good. The agreement for circular, 
inclined orbits is excellent, but this is unsurprising since we used Teukolsky data to fit the circular part of the hybrid 
fluxes. The hybrid scheme also performs very well for orbits of low eccentricity, but the performance is a little worse 
for orbits of higher eccentricity, particularly in the strong-field. Specifically, we see the decrease in \Lz\ with increasing 
eccentricity for p — 2M and p — 3M that was discussed earlier, and find one prograde equatorial orbit that has Lz > 0. 
This is a consequence of the missing eccentricity-dependent terms in the fluxes. In the future, deriving a fit to the 
missing terms using Teukolsky-based data should correct these remaining problems. 
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TABLE I: Comparison of hybrid and Teukolsky-based results. We tabulate the fluxes of energy and angular momentum and 
the rate of change of the inclination angle, t, for a variety of geodesies. These are divided into circular /inclined (top section), 
eccentric/equatorial (middle section) and eccentric/inclined (bottom section). 

C. Sample inspirals 

The main application of the hybrid scheme is to generate inspiral trajectories. These can be used to investigate the 
general properties of extreme mass ratio inspirals, and also as input for the generation of approximate gravitational 
waveforms 26] . Such approximate waveforms are being used for scoping out data analysis issues for the detection of 
EMRIs with LISA 4] and may also find application as detection templates in the final analysis of LISA data. 

In Figure 1111 we illustrate some inspirals in the p — e plane computed using the final form of the hybrid scheme 
(|59|l for a black hole of spin a ^ 0.5M. We show the same inspirals as computed under the original GHK scheme 
for comparison. The improvement of the new scheme over the old one is quite evident and the old pathologies are 
no longer present. In Figure IT^ we show a further sequence of inspirals in the p — e plane computed under the new 
hybrid scheme. These all have initial semi-latus rectum of p — 20M, and have a variety of eccentricities, inclinations 
and spins. Every inspiral has the same basic structure — initially the orbit circularises as it inspirals, but once the 
orbit gets close to plunge, a critical point is reached where e = and the eccentricity then begins to increase until 
the object plunges into the black hole. Orbits of lower initial eccentricity circularise less rapidly. As the spin of the 
central black hole is increased, the plunge point moves closer to the central black hole, and so does the critical radius. 
Orbits of higher inclination plunge further from the central black hole, so the net effect of spin is reduced, as is clear 
from the l ~ 100° panel in Figure [T^ 

Figure 1131 illustrates the same inspirals, but shows how the inspiral proceeds in the p — t plane. We see that, 
typically, l changes only slightly over the inspiral, but always increases. The magnitude of this change increases as 
the spin is increased, and is greater for orbits nearer to polar (i — 90°) than for nearly equatorial (i = 0) orbits. 
Orbits of lower initial eccentricity also exhibit a lower total change in inclination over the inspiral. For low spins, the 
change in t is less than 1°, but it can be somewhat larger (~ 3°) for near polar inspirals into rapidly rotating black 
holes. These features, and the non-zero, but small, increase in inclination is in good agreement with the results of 
Teukolsky-based calculations [l5| . 
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FIG. 11: Inspirals in the {p, e) plane computed using the revised hybrid scheme (solid curves). The dashed curves represent 
the original GHK results and are included for comparison. Initial values for {p, e, l) are (20M, 0.99, 100°), (14M, 0.99, 100°) and 
(20M, 0.4, 100°). The black hole spin is a = 0.5M. 



VIII. CONSERVATIVE SELF-FORCE CORRECTIONS 



In this paper we have focused on improving the computation of the trajectories that EMRIs follow in phase space, 
specifically the evolution of the three principal constants of the motion, E, Lz and Q. These three constants determine 
the shape of any given geodesic. However, a geodesic also has three additional "positional" constants of the motion, 
which essentially label the position of the particle along the geodesic trajectory at some fiducial time. The self-force 
interaction between the inspiralling body and the central black hole acts to change not only the principal constants of 
the motion (the dissipative self- force), but also the positional constants (the conservative self- force). In this analysis 
we have looked only at the dissipative contribution. However, the influence of the conservative self-force on the orbital 
evolution cannot be ignored, since it could lead to several cycles of phase discrepancy in template waveforms over the 
course of an inspiral. In a recent paper |27l| Pound et al. explored a toy problem in which an orbiting charged particle 
experienced both a dissipative and a conservative electromagnetic self-force. They found that in the weak-field the 
phasing contribution from the conservative self-force was comparable to or larger than the dissipative contribution. 
While that particular analysis was not applicable in the strong-field regime and did not look at the gravitational self- 
force directly, it is likely that the conservative gravitational self-force will have a significant influence on the phasing 
of an inspiral. 

The conservative self-force alters the orbit that an inspiralling body with particular energy etc. follows. This change 
in the orbit in principle will modify the radiated fluxes of E, and Q. However, this is a much lower order effect 
and so can be ignored. The expressions presented here can therefore be accurately used to model EMRI phase space 
trajectories. On the other hand, conservative effects should be included when an inspiral path (i.e., the path through 
space which the orbiting body follows) is computed based on the phase space evolution. This is the next step in 
constructing approximate waveforms (see discussion in next section and Figure One way to add conservative 
effects is to include the three positional constants of the motion when parameterising the orbit and evolve these as 
well as the principal constants along the phase space trajectory. This is discussed in more detail in 26] . Perturbative 
calculations have not yet reached the stage of computing the conservative self- force for strong- field orbits, although 
this should be possible in the near future for Schwarzschild orbits There are post-Newtonian expressions for the 
conservative self- force available in the literature (e.g. ,28]), so it should be possible to put together existing results 
to construct a hybrid scheme for the evolution of the positional orbital constants, along similar lines to the hybrid 
scheme described in this paper for the principal constants. However, this is a complicated procedure and we leave it 
for a future paper. 
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FIG. 12: A set of inspirals in the p — e plane computed using the revised hybrid scheme. For each inspiral we initially set 
p — 20M and t — 30° (top panel) ot l — 100° (bottom panel). The inspirals arc shown for two different initial eccentricities, 
e = 0.4 or e = 0.99, and three different values of the central black hole spin, a = 0.3M (solid curves), a = 0.5M (dotted curves) 
and o = 0.9M (dashed curves). 



IX. CONCLUDING REMARKS 



In this paper we have described a number of ways to alter the GHK hybrid scheme in order to expand its vahdity 
and reliabiUty. The inclusion of all of these corrections leads to inspirals that are physically reasonable and seem 
qualitatively correct throughout parameter space. Moreover, direct comparisons of the hybrid fluxes to Teukolsky- 
based results indicate that the approach performs well. However, at present we can only judge the scheme based 
on physical intuition and lessons learnt from studies of specific families of orbits (equatorial-eccentric and circular- 
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FIG. 13: The radiation reaction induced evolution of t during inspiral for the inspirals shown in Figure fT^ For each inspiral 
we initially set p — 20M and l = 30° (top panel) or t = 100° (bottom panel). We used two different values of the initial 
eccentricity, e — 0.4 (dashed curves) or e = 0.99 (solid curves), and three values for the black hole spin (as labelled). 



inclined) and the limited number of fluxes for generic (inclined-eccentric) orbits that are currently available [l4| . As 
more extensive results for generic orbits are generated in the near future, we will be able to more thoroughly test 
this new hybrid scheme and improve it if necessary. An important point to note is that the new scheme is still 
computationally inexpensive, which is why it is valuable to pursue this approach in conjunction with more accurate 
perturbative calculations. 

There are a number of ways in which the scheme might be further improved in the future. The E, and Q flux 
expressions, (|44f) ~ (l45|l and (|55|l . could be improved by including additional terms. The current expressions do not 
include all eccentricity terms, and we saw in section FVI C 21 that additional terms currently available in the literature 
can be easily included, although these do not appear to improve the performance of the hybrid scheme. However, it 
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should hopefully still be possible to gain improvements by using higher order PN results as these become available. 
The spin-independent terms could also be amended by using the 2PN fluxes derived by Gopakumar and Iyer |29] for 
a binary system with non-spinning members moving in quasi-elliptical orbits. A modification like the latter should 
improve the numerical inspirals for e ~ 1. Another approach would be to construct the Fade approximants of the 
PN fluxes, which are typically more effective than the original Taylor-series fluxes Once generic Teukolsky data 
is available, the scheme could also be enhanced by improving the eccentricity-dependent pieces of the E, Lz and L 
fluxes using fits to the Teukolsky results. As discussed in section IVIIII a significant gap in the present hybrid scheme 
is the omission of conservative self-force corrections. While this omission does not affect the phase space inspiral 
trajectories described here, it might significantly influence the gravitational waveform phasing and can not therefore 
be ignored. Using PN results and data from perturbative self-force calculations, it should be possible to include the 
conservative self-force by computing a "positional" phase space trajectory in conjunction with the "principal" phase 
space trajectory described here. This will be pursued in the future. 

A key point to remember in any future version of the hybrid scheme is that the consistency condition H13|l must be 
satisfied by the fluxes in order to produce reasonable inspirals. To incorporate this correction, we added a term to the 
energy flux that is essentially the difference between the prescribed energy flux and what it should be to be consistent 
with the prescribed angular momentum and inclination fluxes. As the flux expressions are pushed to higher and 
higher order, the size of this correction will naturally diminish and we will ultimately converge to the true adiabatic 
inspiral. 

The conditions \i'A\ and (|22|l also have some relevance for Teukolsky-based calculations for generic orbits. Until 
recently, these had been carried out only for circular orbits of arbitrary inclination ^ind equatorial orbits of 
arbitrary eccentricity In both cases, the additional symmetries allow the Carter constant to be evolved correctly. 
More recently, Teukolsky-based calculations have been carried out for "snapshots" of generic orbits Due to the 
difficulties in computing the evolution of the Carter constant for a generic inspiral, these results made use of the 
constant inclination approximation, Z = 0, to compute Q. It is clear from the analysis here that this will not work for 
nearly circular or nearly polar orbits, since we need the correct i in those limits to ensure conditions H13(l and H22|l are 
satisfied. The constant inclination approach is therefore problematic. However, the new formula for the evolution of 
Q presented by Sago et al. |3| allows computation of Q using the same Teukolsky variables used to evaluate E and 
Lz. Future generic codes should make use of this new formula to avoid the consistency problems described here. 

As mentioned earlier, inspirals generated using this hybrid scheme can be used to investigate the qualitative prop- 
erties of EMRIs, and also to construct approximate, "kludge" , EMRI waveforms for use in LISA data analysis. By 
integrating the Kerr geodesic equations along the trajectory through phase space constructed using the method de- 
scribed here, the path followed by a particle can be computed in Boyer-Lindquist coordinates. Identifying these 
coordinates with spherical polar coordinates and constructing the corresponding fiat space quadrupole moment tensor 
allows the generation of an approximate gravitational waveform. Such waveforms are being used to scope out LISA 
data analysis algorithms 4] . In Figure E| we show short snippets of the kludge gravitational waveform from one of 
the inspirals illustrated in Figures [T^ and [TSl with a = 0.9M and initial parameters p — 20M, e = 0.99 and l — 30°. 
The snippets are of 6000s duration (assuming masses of M = lO^M© and ^ = IM©) and are shown at four points 
along the inspiral trajectory. We include this figure for illustration purposes only. More details on the construction 
of kludge waveforms can be found in |2^ . 
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